fn_panel_stocks_PS <- function(df_input, public_info_vars, controls_levels, controls_logs, unlearnable, quarter_lag){
  df <- df_input
  
  if (unlearnable) {
    indep_var_long  <- c("delta_log_payoff_growth", "forecast_growth", "future_forecast_growth")
    indep_var_short <- c("delta_log_payoff_growth", "forecast_growth")
  } else {
    indep_var_long  <- c("delta_log_payoff_growth", "delta_log_payoff_growth_future")
    indep_var_short <- c("delta_log_payoff_growth")
  }
  
  dep_var <- "delta_log_price"
  
  controls <- c()
  controls_short <- c()
  controls_long <- c()
  
  epsilon_vars <- c()
  
  all_vars <- c(controls_levels, controls_logs, public_info_vars)
  
  for (i in 1:length(all_vars)) {
    if(all_vars[i] %in% controls_logs){
      varname <- paste0("log_", all_vars[i])
      varname_sq <- paste0("log_", all_vars[i], "2")
    }
    else if(all_vars[i] %in% public_info_vars){
      varname <- paste0("delta_", all_vars[i])
      varname_sq <- paste0("delta_", all_vars[i], "2")
    }
    else{
      varname <- all_vars[i]
      varname_sq <- paste0(all_vars[i], "2")
    }
    
    varname_lag <- paste0("lagged_", varname)
    varname_lag_sq <- paste0("lagged_", varname_sq)
    
    if(all_vars[i] %in% controls_logs){
      controls <- c(controls, varname_lag)
      epsilon_vars <- c(epsilon_vars, varname_lag, varname_lag_sq)
    }
    else if(all_vars[i] %in% public_info_vars){
      controls <- c(controls, varname)
      epsilon_vars <- c(epsilon_vars, varname, varname_sq)
    }
    else{
      controls <- c(controls, varname_lag)
      epsilon_vars <- c(epsilon_vars, varname_lag, varname_lag_sq)
    }
    
    if(all_vars[i] %in% controls_logs){
      df <- df %>%
        group_by(permno) %>%
        mutate(!!varname := log(get(all_vars[i]))) %>%
        mutate(!!varname_sq := get(varname)^2) %>%
        mutate(!!varname_lag := dplyr::lag(get(varname), n = quarter_lag, default = NA)) %>%
        mutate(!!varname_lag_sq := dplyr::lag(get(varname_sq), n = quarter_lag, default = NA)) %>%
        ungroup()
    }
    else if(all_vars[i] %in% public_info_vars){
      df <- df %>%
        group_by(permno) %>%
        mutate(!!varname := get(all_vars[i]) - dplyr::lag(get(all_vars[i]), n = quarter_lag, default = NA)) %>%
        mutate(!!varname_sq := get(varname)^2) %>%
        ungroup()
    }
    else{
      df <- df %>%
        group_by(permno) %>%
        mutate(!!varname_sq := get(varname)^2) %>%
        mutate(!!varname_lag := dplyr::lag(get(varname), n = quarter_lag, default = NA)) %>%
        mutate(!!varname_lag_sq := dplyr::lag(get(varname_sq), n = quarter_lag, default = NA)) %>%
        ungroup()
    }
    for(j in 1:length(indep_var_long)){
      indep_var <- indep_var_long[j]
      if(all_vars[i] %in% public_info_vars){
        v <- varname
        v_sq <- varname_sq
      }
      else{
        v <- varname_lag
        v_sq <- varname_lag_sq
      }
      
      varname_long <- paste0(v, "X", indep_var)
      
      controls_long <- c(controls_long, varname_long)
      if(indep_var %in% indep_var_short){
        controls_short <- c(controls_short, varname_long)
      }
      df <- df %>%
        mutate(!!varname_long := get(v) * get(indep_var))
    }
  }
  
  for (i in 1:length(all_vars)) {
    if(all_vars[i] %in% controls_logs){
      varname_outer <- paste0("lagged_log_", all_vars[i])
    }
    else if(all_vars[i] %in% public_info_vars){
      varname_outer <- paste0("delta_", all_vars[i])
    }
    else{
      varname_outer <- paste0("lagged_", all_vars[i])
    }
    for (j in i:length(all_vars)) {
      if(j == i){
        next
      }
      if(all_vars[j] %in% controls_logs){
        varname_inner <- paste0("log_", all_vars[j])
      }
      else if(all_vars[j] %in% public_info_vars){
        varname_inner <- paste0("delta_", all_vars[j])
      }
      else{
        varname_inner <- all_vars[j]
      }
      varname <- paste0(varname_outer, "X", varname_inner)
      epsilon_vars <- c(epsilon_vars, varname)
      df <- df %>%
        group_by(permno) %>%
        mutate(!!varname := get(varname_outer)*get(varname_inner)) %>%
        ungroup()
    }
  }
  
  indep_var_long <- c(indep_var_long, controls, controls_long)
  indep_var_short <- c(indep_var_short, controls, controls_short)
  
  df <- df %>% filter(if_all(c(indep_var_long, dep_var), is.finite))
  
  attach(fn_identifying_stocks_PS(df, dep_var, indep_var_long, indep_var_short), warn.conflicts = F)
  
  recovered <- recovered %>%
    mutate(epsilon_long=log(epsilon_long^2),
           epsilon_short=log(epsilon_short^2))
  
  indep_var <- epsilon_vars
  
  recovered <- recovered %>% filter(if_all(c(indep_var), is.finite))
  
  f_long  <- paste("epsilon_long", paste(indep_var,  collapse = " + "), sep = " ~ ")
  f_short <- paste("epsilon_short", paste(indep_var, collapse = " + "), sep = " ~ ")
  
  # Running regressions
  f_long  <- as.formula(f_long)
  f_short <- as.formula(f_short)
  
  reg_long  <- lm(f_long,  data = recovered)
  reg_short <- lm(f_short, data = recovered)
  
  recovered$fitted_long <- exp(reg_long$fitted.values + 0.5*var(reg_long$residuals, na.rm=T))
  recovered$fitted_short <- exp(reg_short$fitted.values + 0.5*var(reg_short$residuals, na.rm=T))
  
  recovered <- recovered %>% mutate(tau_pi_R=(fitted_short-fitted_long)/fitted_short)
  recovered <- recovered %>% mutate(month_end=month(date_fred),
                                    year_end=year(date_fred))
  
  return(list(results=recovered, reg_long_eps=reg_long, reg_short_eps=reg_short,
              reg_long_PS=reg_long_PS, reg_short_PS=reg_short_PS))
}